#pragma once

#include <cfloat>
#include <cmath>
#include <functional>
#include <initializer_list>
#include <iostream>
#include <iomanip>
#include <string>
#include <vector>
#include <memory>
#include <stdexcept>
#include <cstdlib>
#include <bits/stdc++.h>

using Real = double;

/**
 * @file Tensor.h
 * @author xingyifan
 * @date 2022-07-30 19:26
 *
 * @description: 整型算法
 */

namespace IntOp
{
    // 指标加法
    template <int Dim>
    inline std::shared_ptr<int> add(const int *ind1, const int *ind2)
    {
        // 创建临时指标
        std::shared_ptr<int> tindex = std::shared_ptr<int>(new int[Dim]);
        for (int i = 0; i < Dim; i++)
        {
            tindex.get()[i] = ind1[i] + ind2[i];
        }
        return tindex;
    }

    // 整型指数函数
    inline int power(int n, int d)
    {
        int s = 1;
        for (int i = 0; i < d; i++)
        {
            s *= n;
        }
        return s;
    }

    // 比较函数
    inline int min(int a, int b)
    {
        return (a < b) ? a : b;
    }

    // 比较函数
    inline int max(int a, int b)
    {
        return (a > b) ? a : b;
    }

    // 绝对值函数
    inline int abs(int a)
    {
        return (a > 0) ? a : -a;
    }
};

/**
 * @file Tensor.h
 * @author xingyifan
 * @date 2022-07-29 18:28
 *
 * @description: 检测异常
 */

class ErrorTest
{
public:
    /**
     * @file Tensor.h
     * @author xingyifan
     * @date 2022-07-29 18:41
     *
     * @description: 根据表达式是否错误判断是否抛出异常
     */

    static void Test(bool expr, int lineNum, std::string file, std::string func, std::string warning)
    {
        try
        {
            // expr 错误则抛出异常
            if (!expr)
                throw ErrorTest(lineNum, file, func, warning);
        }
        catch (const ErrorTest &e)
        {
            // 输出异常信息，暂停程序
            e.what();
            std::system("pause");
        }
    }

private:
    ErrorTest(int lineNum, std::string file, std::string func, std::string warning)
        : lineNum_(lineNum), file_(file), func_(func), warning_(warning) {}

    /**
     * @file Tensor.h
     * @author xingyifan
     * @date 2022-07-29 18:38
     *
     * @description: 输出错误信息
     */

    void what() const
    {
        std::cout << file_ << ":" << lineNum_ << ": error: In function \"" << func_ << "\", " << warning_ << std::endl;
        std::cout << std::endl;
    }

private:
    std::string warning_; // 错误提示
    std::string file_;    // 错误文件名
    std::string func_;    // 错误函数名
    int lineNum_;         // 错误行号
};

/**
 * @file Tensor.h
 * @author xingyifan
 * @date 2022-07-25 12:59
 *
 * @description: 张量类型，注意使用 Tensor 需要时刻关注指标的刷新问题！！！！！
 */

template <int Dim>
class Tensor
{
public:
    /**
     * @file Tensor.h
     * @author xingyifan
     * @date 2022-07-25 12:54
     *
     * @description: 构造函数
     */

    // 初始宽度为 1
    Tensor(int base = 1) : x_(new Real[IntOp::power(base, Dim)]), size_(IntOp::power(base, Dim)), base_(base)
    {
        // 初始化所有索引
        for (int i = 0; i < Dim; ++i)
        {
            index_[i] = 0;
        }
    }
    // 基数+填充值
    Tensor(int base, Real r) : Tensor(base)
    {
        int n = size();
        for (int i = 0; i < n; ++i)
        {
            x_.get()[i] = r;
        }
    }
    // 基数+初始化列表
    Tensor(int base, std::initializer_list<Real> l) : Tensor(base)
    {
        int i = 0;
        for (auto it = l.begin(); it != l.end(); i++)
        {
            x_.get()[i] = *it++;
        }
    }
    // 初始化列表，方便向量初始化
    Tensor(std::initializer_list<Real> l) : Tensor(l.size())
    {
        int i = 0;
        for (auto it = l.begin(); it != l.end(); i++)
        {
            x_.get()[i] = *it++;
        }
    }
    // 拷贝构造函数
    Tensor(const Tensor &t) : Tensor(t.width())
    {
        int n = size();
        for (int i = 0; i < n; i++)
        {
            x_.get()[i] = t[i];
        }
    }
    // 移动构造函数
    Tensor(const Tensor &&t) : x_(std::move(t.x_)), size_(t.size()), base_(t.base())
    {
        for (int i = 0; i < Dim; ++i)
        {
            index_[i] = 0;
        }
    }
    // 不同阶张量移动构造
    template <int D>
    Tensor(int base, const Tensor<D> &&t) : x_(std::move(t.smartptr())), size_(t.size()), base_(base)
    {
        for (int i = 0; i < Dim; ++i)
        {
            index_[i] = 0;
        }
    }

    /**
     * @file Tensor.h
     * @author xingyifan
     * @date 2022-07-25 12:54
     *
     * @description: 赋值语句
     */

    // 拷贝赋值语句
    void operator=(const Tensor &t)
    {
        // 改变基和容器大小
        base_ = t.base_;
        size_ = t.size();

        // 创建新容器
        x_ = std::shared_ptr<Real>(new Real[size_]);

        // 赋值
        for (int i = 0; i < size_; i++)
        {
            x_.get()[i] = t[i];
        }
    }
    // 移动赋值语句
    void operator=(const Tensor &&t)
    {
        size_ = t.size();
        base_ = t.base_;
        x_ = std::move(t.x_);
    }
    // 允许初始化列表赋值
    void operator=(std::initializer_list<Real> l)
    {
        // 重新申请内存
        size_ = l.size();
        x_ = std::shared_ptr<Real>(new Real[size_]);

        // 赋值
        int i = 0;
        for (auto it = l.begin(); it != l.end(); i++)
        {
            x_.get()[i] = *it++;
        }
    }
    // 允许数组赋值
    void operator=(Real *l)
    {
        int n = size();
        for (int i = 0; i < n; i++)
        {
            x_.get()[i] = l[i];
        }
    }
    // 不同阶移动赋值
    template <int D>
    void move(int base, const Tensor<D> &&t)
    {
        size_ = t.size();
        base_ = base;
        x_ = std::move(t.smartptr());
    }

    /**
     * @file Tensor.h
     * @author xingyifan
     * @date 2022-07-25 12:54
     *
     * @description: 指标操作
     */

    // 标准移动指标（包含进位操作）
    bool operator++(int)
    {
        int ind = 0;

        while (index_[ind] == base_ - 1 && ind + 1 < Dim)
        {
            index_[ind] = 0;
            ind++;
        }
        index_[ind] = (index_[ind] + 1) % base_;

        // 返回 false 说明到达边界
        return (ind != Dim - 1 || index_[ind] != 0);
    }

    // 固定指定维的指标，在其它方向上移动，对于一阶张量无效
    bool move(int d)
    {
        int ind = 0;
        int ind2 = (ind >= d) ? ind + 1 : ind;

        // 如果达到 base_ - 1，此位设为零
        while (ind + 1 < Dim - 1 && index_[ind2] == base_ - 1)
        {
            // 进位，此位归零
            index_[ind2] = 0;
            ind++;

            // 移动标记
            ind2 = (ind >= d) ? ind + 1 : ind;
        }

        // 进位加一
        index_[ind2] = (index_[ind2] + 1) % base_;

        // 调整边界
        int D = (d == Dim - 1) ? d - 1 : Dim - 1;

        // 返回 false 说明到达边界
        return (ind2 != D || index_[ind2] != 0);
    }

    // 在指定维上移动，可循环
    void forward(int d) { index_[d] = (index_[d] + 1) % base_; }
    void backward(int d) { index_[d] = (index_[d] - 1 + base_) % base_; }

    // 刷新指标
    void fresh()
    {
        for (int i = 0; i < Dim; i++)
        {
            index_[i] = 0;
        }
    }

    // 得到细网格上对应的指标
    void fine(Tensor &fineTensor)
    {
        for (int d = 0; d < Dim; d++)
        {
            fineTensor.index_[d] = index_[d] * 2;
        }
    }

    // 传递指标
    void passIndex(Tensor &t)
    {
        for (int d = 0; d < Dim; d++)
        {
            t.index_[d] = index_[d];
        }
    }

    // 传入指标，计算总指标
    int totalIndex(int *index = nullptr)
    {
        if (index == nullptr)
            index = index_;
        int ind = 0;
        for (int d = 0; d < Dim; d++)
        {
            ind += index[d] * IntOp::power(base_, d);
        }
        return ind;
    }

    /**
     * @file Tensor.h
     * @author xingyifan
     * @date 2022-07-25 12:55
     *
     * @description: 获得成员
     */

    int size() const { return size_; }  // 容器大小 = 宽度^维数
    int width() const { return base_; } // 返回容器宽度
    int dim() const { return Dim; }     // 返回张量阶

    // 返回 d 维分量指标
    int index(int d) const
    {
        // 数组越界预警
        ErrorTest::Test(d >= 0 && d < Dim, __LINE__, __FILE__, __func__, "You are out of range of the Index!");
        return index_[d];
    }
    // 允许修改 d 维分量指标
    int &index(int d)
    {
        // 数组越界预警
        ErrorTest::Test(d >= 0 && d < Dim, __LINE__, __FILE__, __func__, "You are out of range of the Index!");
        return index_[d];
    }

    // 允许修改指标基
    int base() const { return base_; } // 返回指标基
    int &base() { return base_; }      // 返回指标基的引用

    // 获得头部智能指针
    const std::shared_ptr<Real> &smartptr() const { return x_; }

    /**
     * @file Tensor.h
     * @author xingyifan
     * @date 2022-07-28 23:39
     *
     * @description: 设置预警，越界立即报错并暂停程序
     */

    // 通过一维指标访问和修改元素
    Real operator[](int i) const
    {
        // 数组越界预警
        ErrorTest::Test(i >= 0 && i < size(), __LINE__, __FILE__, __func__, "You are out of range of the Tensor!");
        return x_.get()[i];
    }
    Real &operator[](int i)
    {
        // 数组越界预警
        ErrorTest::Test(i >= 0 && i < size(), __LINE__, __FILE__, __func__, "You are out of range of the Tensor!");
        return x_.get()[i];
    }

    // 通过多重指标访问元素
    Real value(int step = 0) const
    {
        int i = totalIndex() + step;

        // 数组越界预警
        ErrorTest::Test(i >= 0 && i < size(), __LINE__, __FILE__, __func__, "You are out of range of the Tensor!");
        return x_.get()[i];
    }
    Real &value(int step = 0)
    {
        int i = totalIndex() + step;

        // 数组越界预警
        ErrorTest::Test(i >= 0 && i < size(), __LINE__, __FILE__, __func__, "You are out of range of the Tensor!");
        return x_.get()[i];
    }

    // 获得当前元素指针
    Real *ptr()
    {
        int i = totalIndex();

        // 数组越界预警
        ErrorTest::Test(i >= 0 && i < size(), __LINE__, __FILE__, __func__, "You are out of range of the Tensor!");
        return &(x_.get()[i]);
    }

    /**
     * @file Tensor.h
     * @author xingyifan
     * @date 2022-07-25 12:55
     *
     * @description: 类型转换
     */

    // 显式转换为实数，使用隐式转换容易导致歧义
    explicit operator Real() const { return x_.get()[0]; }

    // 隐式转换为向量
    operator std::vector<Real>() const
    {
        int n = size();
        std::vector<Real> res(n);
        for (int i = 0; i < n; ++i)
        {
            res[i] = x_.get()[i];
        }
        return res;
    }

    /**
     * @file Tensor.h
     * @author xingyifan
     * @date 2022-07-25 12:56
     *
     * @description: 获得不同阶的 Tensor
     */

    // 特殊的张量积，t1 作为新的分量补充到当前张量
    template <int D>
    Tensor<Dim + D> product(const Tensor<D> &t1)
    {
        Tensor<Dim + D> t(base_);
        int n = size();
        int n1 = t1.size();
        int step = IntOp::power(base_, Dim);

        // 循环增加指标
        for (int i = 0; i < n1; ++i)
        {
            for (int j = 0; j < n; ++j)
            {
                t[i * step + j] = x_.get()[j] * t1[i];
            }
        }
        return t;
    }

    // 返回 D 阶 Tensor
    template <int D>
    Tensor<D> component(int i) const
    {
        // 数组越界预警
        ErrorTest::Test(IntOp::power(base_, D) * i < size_, __LINE__, __FILE__, __func__, "You are out of range of the Tensor!");

        // 步长
        int step = IntOp::power(base_, Dim - D);

        Tensor<D> t(base_);
        int n = t.size();
        for (int j = 0; j < n; j++)
        {
            t[j] = *(x_.get() + step * j + i);
        }
        return t;
    }

    // 获得包围当前值的所有值
    Tensor<Dim> cell()
    {
        // 模拟指标，每个分量 -1
        int tindex[Dim];
        for (int i = 0; i < Dim; i++)
        {
            tindex[i] = index_[i] - 1;
        }

        // 检测对应点是否在边界上，只要在边界上就报错
        for (int d = 0; d < Dim; d++)
        {
            if (tindex[d] == 0 || tindex[d] == base_ - 1)
            {
                ErrorTest::Test(false, __LINE__, __FILE__, __func__, "You are out of range of the Tensor!");
                break;
            }
        }

        // 返回容器
        Tensor<Dim> val(3);

        // 计算并填入指标对应值
        do
        {
            // 获得加和指标
            std::shared_ptr<int> index = IntOp::add<Dim>(tindex, val.index_);
            val.value() = *(x_.get() + totalIndex(index.get()));
        } while (val++);

        return val;
    }

    /**
     * @file Tensor.h
     * @author xingyifan
     * @date 2022-07-25 12:57
     *
     * @description: 边界判定
     */

    // 判断此方向上是上边界 1 还是下边界 -1 (均为 true)，或者不是边界 0 (false)
    int isBound(int d) const
    {
        if (index_[d] == 0)
            return -1;
        if (index_[d] == base_ - 1)
            return 1;
        return 0;
    }

    // 判断此方向上是靠近上边界 2,1 还是靠近下边界 -2,-1 (均为 true)，或者不靠近边界 0 (false)
    int nearBound(int d) const
    {
        if (index_[d] == 0)
            return -2;
        if (index_[d] == 1)
            return -1;

        if (index_[d] == base_ - 2)
            return 1;
        if (index_[d] == base_ - 1)
            return 2;
        return 0;
    }

    // 查找边界分量，不在边界返回 -1
    int findBound() const
    {
        for (int d = 0; d < Dim; d++)
        {
            if (isBound(d))
                return d;
        }
        return -1;
    }

    /**
     * @file Tensor.h
     * @author xingyifan
     * @date 2022-07-25 12:57
     *
     * @description: 张量四则远算，需要注意添加 const 修饰，否则 const Tensor 不能调用这些函数
     */

// 关于数的四则运算 + - * /
#define TENSOR_BINARY_OP(OpName, Op) \
    Tensor OpName(Real r) const      \
    {                                \
        Tensor t = *this;            \
        int n = size();              \
        for (int i = 0; i < n; i++)  \
        {                            \
            t[i] = t[i] Op r;        \
        }                            \
        return t;                    \
    }

    TENSOR_BINARY_OP(operator+, +)
    TENSOR_BINARY_OP(operator-, -)
    TENSOR_BINARY_OP(operator*, *)
    TENSOR_BINARY_OP(operator/, /)
#undef TENSOR_BINARY_OP

// 关于数的四则运算 += -= *= /=
#define TENSOR_BINARY_OP(OpName, Op) \
    void OpName(Real r) const        \
    {                                \
        int n = size();              \
        for (int i = 0; i < n; i++)  \
        {                            \
            x_.get()[i] Op r;        \
        }                            \
    }

    TENSOR_BINARY_OP(operator+=, +=)
    TENSOR_BINARY_OP(operator-=, -=)
    TENSOR_BINARY_OP(operator*=, *=)
    TENSOR_BINARY_OP(operator/=, /=)
#undef TENSOR_BINARY_OP

// 关于张量的四则运算 + - * /
#define TENSOR_BINARY_OP(OpName, Op, Ref)    \
    Tensor OpName(const Tensor Ref t0) const \
    {                                        \
        Tensor t = *this;                    \
        int n = size();                      \
        for (int i = 0; i < n; i++)          \
        {                                    \
            t[i] = t[i] Op t0[i];            \
        }                                    \
        return t;                            \
    }

    TENSOR_BINARY_OP(operator*, *, &)
    TENSOR_BINARY_OP(operator/, /, &)
    TENSOR_BINARY_OP(operator+, +, &)
    TENSOR_BINARY_OP(operator-, -, &)
#undef TENSOR_BINARY_OP

// 关于张量的四则运算 += -= *= /=
#define TENSOR_BINARY_OP(OpName, Op, Ref)  \
    void OpName(const Tensor Ref t0) const \
    {                                      \
        int n = size();                    \
        for (int i = 0; i < n; i++)        \
        {                                  \
            x_.get()[i] Op t0[i];          \
        }                                  \
    }

    TENSOR_BINARY_OP(operator*=, *=, &)
    TENSOR_BINARY_OP(operator/=, /=, &)
    TENSOR_BINARY_OP(operator+=, +=, &)
    TENSOR_BINARY_OP(operator-=, -=, &)
#undef TENSOR_BINARY_OP

    // 取负值
    friend Tensor operator-(Tensor &t) { return t * -1; }

    // 张量乘法：Dim - D 阶张量 = Dim 阶张量 * D 阶张量
    template <int D>
    Tensor<Dim - D> operator*(const Tensor<D> &t2) const
    {
        Tensor<Dim - D> t(base_);
        int n = t2.size();
        for (int i = 0; i < n; ++i)
        {
            t = t + component<Dim - D>(i) * t2[i];
        }
        return t;
    }

    /**
     * @file Tensor.h
     * @author xingyifan
     * @date 2022-07-25 12:57
     *
     * @description: 输出函数
     */

    friend std::ostream &operator<<(std::ostream &out, const Tensor &t)
    {
        int ind = 0;
        int n = t.size();
        while (n > ind)
        {
            for (int i = 0; i < t.base_; i++, ind++)
            {
                out << std::setw(5) << t[ind] << " ";
            }
            out << std::endl;
        }
        return out;
    }

private:
    std::shared_ptr<Real> x_; // 容器
    int size_;                // 容器大小
    int base_;                // 指标基，即容器宽度
    int index_[Dim];          // 多重指标
};

/**
 * @file Tensor.h
 * @author xingyifan
 * @date 2022-07-30 19:24
 *
 * @description: 特殊算法
 */

namespace TensorOp
{
    // 特殊张量指数
    template <int Dim>
    inline Tensor<Dim> power(const Tensor<1> &t)
    {
        int n = t.size();

        // 构造两个容器
        Tensor<Dim> t1(n);
        Tensor<Dim> t2(n);

        // 初始化其中一个
        for (int i = 0; i < n; i++)
        {
            t1[i] = t[i];
        }

        // 循环乘积
        for (int d = 1; d < Dim; d++)
        {
            int step = IntOp::power(n, d);
            for (int i = 0; i < step; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    t2[i * n + j] = t1[i] * t[j];
                }
            }
            for (int i = 0; i < step * n; i++)
            {
                t1[i] = t2[i];
            }
        }
        return t1;
    }

    // 张量点积
    template <int Dim>
    inline Real dot(const Tensor<Dim> &t1, const Tensor<Dim> &t2)
    {
        int n = t1.size();
        Real sum = 0;
        for (int i = 0; i < n; ++i)
        {
            sum += t1[i] * t2[i];
        }
        return sum;
    }

    // 张量范数：直接使用向量范数
    template <int Dim>
    inline Real norm(const Tensor<Dim> &t, int p = 0)
    {
        // 张量规模
        int n = t.size();

        // 存放结果
        Real s = 0;
        if (p == 2)
        {
            for (int i = 0; i < n; i++)
            {
                s += t[i] * t[i];
            }
            s = sqrt(s);
        }
        else if (p == 1)
        {
            for (int i = 0; i < n; i++)
            {
                s += fabs(t[i]);
            }
        }
        else
        {
            for (int i = 0; i < n; i++)
            {
                Real d = fabs(t[i]);
                if (s < d)
                {
                    s = d;
                }
            }
        }
        return s;
    }

    // 使用 Grid function 计算范数
    template <int Dim>
    inline Real grid(const Tensor<Dim> &t, int p = 0)
    {
        // 网格规模
        int n = t.size();

        // 分割宽度
        Real h = 1.0 / (t.width() - 1);

        // 存放结果
        Real s = 0;
        if (p == 2)
        {
            for (int i = 0; i < n; i++)
            {
                s += t[i] * t[i];
            }
            s = sqrt(pow(h, Dim) * s);
        }
        else if (p == 1)
        {
            for (int i = 0; i < n; i++)
            {
                s += pow(h, Dim) * fabs(t[i]);
            }
        }
        else
        {
            for (int i = 0; i < n; i++)
            {
                Real d = fabs(t[i]);
                if (s < d)
                {
                    s = d;
                }
            }
        }
        return s;
    }

    /**
     * @file Tensor.h
     * @author xingyifan
     * @date 2022-08-30 16:06
     *
     * @description: 利用移动构造优化运算
     */

#define TENSOR_BINARY_OP(OpName, Op, Ref)                                \
    template <int Dim>                                                   \
    inline Tensor<Dim> OpName(Tensor<Dim> Ref t1, const Tensor<Dim> &t2) \
    {                                                                    \
        Tensor<Dim> t = std::move(t1);                                   \
        int n = t.size();                                                \
        for (int i = 0; i < n; i++)                                      \
        {                                                                \
            t[i] = t[i] Op t2[i];                                        \
        }                                                                \
        return t;                                                        \
    }

    TENSOR_BINARY_OP(add, +, &)
    TENSOR_BINARY_OP(add, +, &&)
    TENSOR_BINARY_OP(sub, -, &)
    TENSOR_BINARY_OP(sub, -, &&)
#undef TENSOR_BINARY_OP

#define TENSOR_BINARY_OP(OpName, Op, Ref)                 \
    template <int Dim>                                    \
    inline Tensor<Dim> OpName(Tensor<Dim> Ref t1, Real r) \
    {                                                     \
        Tensor<Dim> t = std::move(t1);                    \
        int n = t.size();                                 \
        for (int i = 0; i < n; i++)                       \
        {                                                 \
            t[i] = t[i] Op r;                             \
        }                                                 \
        return t;                                         \
    }

    TENSOR_BINARY_OP(multiple, *, &)
    TENSOR_BINARY_OP(multiple, *, &&)
    TENSOR_BINARY_OP(divide, /, &)
    TENSOR_BINARY_OP(divide, /, &&)
#undef TENSOR_BINARY_OP

    /**
     * @file Tensor.h
     * @author xingyifan
     * @date 2022-07-25 12:58
     *
     * @description: 收敛阶计算
     */

    template <int Dim>
    inline Real rate(const Tensor<Dim> &t)
    {
        int n = t.size();
        Real r = 0;
        for (int i = 0; i < n - 1; i++)
        {
            r += log2(t[i] / t[i + 1]) / (n - 1);
        }
        return r;
    }

    // 对数张量
    template <int Dim>
    inline Tensor<Dim> log(const Tensor<Dim> &t)
    {
        int n = t.size();
        Tensor<Dim> logt(t.base());
        for (int i = 0; i < n; i++)
        {
            logt[i] = log2(t[i]);
        }
        return logt;
    }
};

/**
 * @file Tensor.h
 * @author xingyifan
 * @date 2022-07-26 18:26
 *
 * @description: 标准容器
 */

using Vector = Tensor<1>;
using Matrix = Tensor<2>;

/**
 * @file Tensor.h
 * @author xingyifan
 * @date 2022-07-26 18:27
 *
 * @description: 函数包装
 */

using Func = std::function<Vector(Real)>;           // 一元向量函数
using FuncX = std::function<Vector(Vector)>;        // 多元向量函数
using FuncXt = std::function<Vector(Vector, Real)>; // 多元 + t 向量函数